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Abstract 

Pariser-Parr-Pople (P-P-P) model Hamiltonian is employed frequently to study 
the electronic structure and optical properties of 7r-conjugated systems. In this 
paper we describe a Fortran 90 computer program which uses the P-P-P model 
Hamiltonian to solve the Hartree-Fock (HF) equation for infinitely long, one- 
dimensional, periodic, 7r-electron systems. The code is capable of computing 
the band structure, as also the linear optical absorption spectrum, by using 
the tight-binding (TB) and the HF methods. Furthermore, using our program 
the user can solve the HF equation in the presence of a finite external electric 
field, thereby, allowing the simulation of gated systems. We apply our code 
to compute various properties of polymers such as irans-polyacetylene (i-PA), 
poly-para-phenylene (PPP), and armchair and zigzag graphene nanoribbons, in 
the infinite length limit. 
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of 64-bit Fedora including Fedora 14 (kernel version 2.6.35.12-90) 
Programming language used: Fortran 90 

Compilers used: Program has been tested with Intel Fortran Compiler (non- 
commercial version 11.0.074) and gfortran compiler (gcc version 4.5.1) with 
optimization option -O. 

Libraries needed: This program needs to link with LAPACK/BLAS libraries 
compiled with the same compiler as the program. For the Intel Fortran Com- 
piler we used the ACML library version 4.4.0, while for the gfortran compiler 
we used the libraries supplied with the Fedora distribution. 
Number of bytes in distributed program, including test data, etc.: ... size of the 
tar file 

Number of lines in distributed program, including test data, etc.: ... lines in the 
tar file 

Card punching code: ASCII 

Nature of physical problem: The electronic structure of one-dimensional peri- 
odic TT-conjugated systems is an intense area of research at present because of 
the tremendous interest in the physics of conjugated polymers and graphene 
nanoribbons. The computer program described in this paper provides an ef- 
ficient way of solving the Hartree-Fock equations for such systems within the 
P-P-P model. In addition to the Bloch orbitals, band structure, and the density 
of states, the program can also compute quantities such as the linear absorption 
spectrum, and the electro-absorption spectrum of these systems. 
Method of Solution: For a one-dimensional periodic 7r-conjugated system ly- 
ing in the xy-plane, the single-particle Bloch orbitals are expressed as linear 
combinations of p^-orbitals of individual atoms. Then using various parameters 
defining the P-P-P Hamiltonian, the Hartree-Fock equations are set up as a ma- 
trix eigenvalue problem in the fc-space. Thereby, its solutions are obtained in a 
self-consistent manner, using the iterative diagonalizing technique at several fc- 
points. The band structure and the corresponding Bloch orbitals thus obtained 
are used to perform a variety of calculations such as the density of states, linear 
optical absorption spectrum, electro-absorption spectrum, etc. 
Running Time: Most of the examples provided take only a few seconds to run. 
For a large system, however, depending on the system size, the run time may 
be a few minutes to a few hours. 
Unusual features of the program: None 



1. Introduction 

Conjugated molecules and polymers have been actively investigated theoret- 
ically as well as experimentally in recent years, (author?) [l|, [3,13 because of 
their potential applications in manufacture of optoelectronic devices, and solar 
cells (author?) |4j]. The low- lying excitations in such materials are characterized 
by the tt electrons, which have itinerant nature, and form the energy levels near 
the chemical potential (Fermi level). Recently, the field of 7r-conjugated systems 
has received a tremendous boost with the synthesis of graphene(author?) 
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and its heterostructures such as graphene nanoribbons(author?) which ex- 
hibit exotic transport and electronic properties, leading to the possibility of 
future electronic devices based upon graphene rather than silicon(author?) 
0, S, 01- Because of these recent advances, theoretical studies of 7r-electron 
systems have come to the forefront of physics. Most of the theoretical meth- 
ods used for describing the electronic structure of these materials can be clas- 
sified as: (a) fully ab initio approaches based upon the mean-field methods 
such as the density-functional theory (DFT) (author?) Ep, El; E H Q 
[lit or the Hartree-Fock (HF) method(author?) [H, [ll [ij and (b) ap- 
proaches based upon effective 7r-electron models such as the tight-binding (TB) 



model(author?) |20|, l2l|, l22| , and the Hubbard(author?) [23|, l2J, l25| or the 
extended Hubbard mo del (author?) [2^. The ab initio methods are generally 
computationally intensive because they make no distinction between the a and 
the TT electrons of the system, and therefore require the use of large basis sets to 
provide a reasonable description of the electronic structure of such systems. In 
case of graphene nanoribbons (GNRs) of large widths, and also polymers with 
large unit cells, the number of degrees of freedom involved in the problem may 
impose severe limitations on the problems which can be solved computationally. 
On the other hand the main advantage of the effective 7r-electron model Hamil- 
tonians is that they explicitly deal only with the tt electrons, thereby reducing 
the degrees of freedom considerably, and, thus allowing the simulation of much 
larger systems as compared to the ab initio approaches. Their disadvantage, of 
course, is that they are semiempirical in nature, and, therefore, the parameters 
involved in them are often determined using the spectroscopic data of a suitable 
model system. Nevertheless, when calculations on very large systems need to 
be performed, often it is virtually impossible to use the ab initio approaches, 
and, therefore, model Hamiltonians provide an attractive alternative. Even 
for smaller systems, model Hamiltonians allow us to understand the underly- 
ing physics in simple terms, therefore, aforesaid 7r-electron approaches are very 
popular when it comes to calculations of the electronic structure of graphene 
and its nanostructures. 

Among the effective 7r-electron approaches, the TB model (called the Hiickel 
model in the chemistry literature) is the simplest, but it does not incorporate the 
effect of electron-electron (e-e) repulsion. The Hubbard model and its extended 
versions go beyond the TB approach, and incorporate the short-range parts of 
the Coulomb repulsion such as the on-site, and the nearest-neighbor interac- 
tions, respectively. In chemistry literature it is well-known that in 7r-electron 
systems such as various aromatic molecules and conjugated polymers, the role 
of long-range e-e interactions cannot be ignored when it comes to their elec- 
tronic structure(author?) [H, 13, 3 ■ And, indeed, Pariser-Parr-Pople (P-P-P) 
model Hamiltonian(author?) [27|, which is an effective 7r-electron Hamiltonian 
incorporating long-range e-e interaction, has been used with considerable suc- 
cess in describing the physics of such systems (author?) P, 0, Thus, it is 
logical to conclude that such e-e interactions will also be important in under- 
standing the physics of the graphene- based materials. In recent years, o ur group 
and collaborators (author?) ^ H M, M, M, S S, M, M, 113, H M, 
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along with numerous other groups(author?) P, [3, H, \M, \M, S, H H Sll, 



have used the P-P-P model to study the electronic structure and optical prop- 
erties of conjugated molecules and oligomers. For finite 7r-conjugated systems, 
in our group, we have developed a general-purpose HF program employing the 
P-P-P model, which is available to anyone for scientific work(author?) (47| . 
However, GNRs and some 7r-conjugated polymers are believed to be quasi-one- 
dimensional systems, which need to be studied in the infinite-length limit. Thus, 
in order to study their electronic structure and related properties, using the P- 
P-P model (or any other Hamiltonian) , one needs to impose periodic boundary 
conditions, and perform infinite lattice sums to account for their infinite extent, 
which, our earlier computer program(author?) lacks. With this aim in 
mind, we recently developed a P-P-P model-based computer program, which 
solves both the restricted HF (RHF) and the unrestricted HF (UHF) equa- 
tions for one-dimensional periodic systems, and used it to study the electronic 
structure and optical properties of mono-layer and multilayer GNRs of various 
kinds(author?) (isl. H^. The ability to solve the UHF equations allows us to 
explore magnetic properties of polymers and GNRs. The aim of the present pa- 
per is to describe the computer program in detail, and make it available for use 
by anyone interested in the physics of GNRs and 7r-conjugated polymers. The 
program is capable of computing the total energy, the band structure, the den- 
sity of states (DOS), and also the interband optical absorption spectrum in form 
of the frequency-depdendent dielectric response tensor, both with and without 
an external homogeneous electric field. The fact that our program can solve 
the HF equations in the presence of an external electric field allows the user to 
explore gated configurations of polymers and GNRs. In this work, we demon- 
strate the capabilities of our program by performing calculations on polymers 
<rans-polyacetylene (t-PA), poly-para-phenylene (PPP), armchair GNRs (AG- 
NRs), and zigzag GNRs (ZGNRs) of various widths. We also note that Nakada 
et aZ. (author?) js^ reported a P-P-P model based band structure calculation 
of both AGNRs and ZGNRs several years back. 

The remainder of the paper is organized as follows. In section [5] we briefly 
review the theory associated with the P-P-P model Hamiltonian. Next, in 
section [3] we discuss the general structure of our computer program, and also 
describe its constituent subroutines. In section |4] we briefly describe how to 
install the program on a given computer system, and to prepare the input files. 
Results of various example calculations using our program are presented and 
discussed in section [5l Finally, in section IH we present our conclusions, as well 
as discuss possible future directions. 

2. Theory 

In this section we briefiy discuss the theory behind our computer program. 
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2.1. Pariser-Parr-Pople Hamiltonian 

The P-P-P model Hamiltonian(author?) [l^l, with one 7r-electron per car- 
bon atom, is given by 

U Ui^n.i + Vij (n, - l){nj - 1) (1) 

i i<j 

where represents the site energy associated with the ith carbon atom, 
creates an electron of spin a on the orbital of atom i, ni^ = c\^Cicj is the 
number of electrons with the spin a, and Ui = niu is the total number of 
electrons on atom i. The parameters U and Vij are the on-site and long-range 
Coulomb interactions, respectively, while tij is the one-electron hopping matrix 
element. On setting Vij = (with U ^ Q), the Hamiltonian reduces to the 
Hubbard model, while on setting both ?7 = and Vij = 0, the tight-binding 
(TB) model is obtained. The parametrization of Coulomb interactions is Ohno 
like (author?) [5l|, 

" «,,,(! + 0.6117i?2^.)'/2 ' 

where, Kij depicts the dielectric constant of the system which can simulate the 
effects of screening, and Ri,j is the distance in A between the i-th and the 
j-th carbon atoms. In our earlier work on GNRs(author?) we used the 
ab-initio GW band structure of mono layer AGNR-12 (AGNR-A'yi, denotes an 
AGNR with Na dimer lines across the width) reported by Son et aZ. (author?) 

to obtain a set of "modified screened Coulomb parameters," with U = 6.0 
eV and j = 2.0 (i ^ j) and Ki^i = 1, which are slightly different from the 
screened parameters reported initially by Chandross and Mazumdar (author?) 



52| . with U ~ 8.0 eV and Kij — 2.0 (i ^ j) and k.^ j = 1, aimed at describing 
the optical properties of phenyl-based polymers within the P-P-P model. The 
modified screened parameters provided good agreement between our HE band 
gaps, and those obtained by the GW method for a few AGNRs, however, for 
ZGNRs the agreement was not good(author?) [48|. In this work, we examine 
the issue of the choice of Coulomb parameters in a critical manner, and conclude 
that no single set of parameters gives uniformly good agreements between our 
results and the GW results for all types of GNRs. In section \5\ where this 
issue is investigated, we find that a set of parameters which provides excellent 
agreement between ours and GW results for a class of GNRs, may simply fail 
to reproduce such agreement for another class of GNRs. In other words, the 
choice of Coulomb parameters which will lead to good agreement between our 
HP results, and the ab initio GW ones, depends upon the geometries of the 
GNRs in question. Thus, after trying a number of Coulomb parameters, and 
in the absence of any experimental data on the band gaps of GNRs, we have 
decided it is best to use the original screened parameters of Chandross and 
Mazumdar (author?) [H^l, with U = 8.0 eV and Kij = 2.0 {i ^ j) and KiA = 1 
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for all the GNRs. But, we would like to emphasize that the user has the freedom 
to choose a different set of Coulomb parameters as per the requirements, and, 
we encourage such experimentation. 

2.2. Unrestricted Hartree-Fock Equations 

We have implemented the RHF and the UHF methods within the P-P- 
P model, using the standard linear combination of atomic orbitals (LCAO) 
formalism. We shall review the basics of the formalism for the UHF method 
(also known as the Pople-Nesbet formalism(author?) [H^), from which the 
corresponding equations for the RHF method can be easily deduced. In this 
approach, the n— th Bloch orbital of the system corresponding to the spin up 
electrons (a//3 will denote spin up/down electrons) is expressed as a linear 
combination of m basis functions per unit cell, 

rn 

^i"Hk) = Y.^lnik)Mk) (3) 

where C^n^(fc)'s represent the linear expansion coefficients, to be determined 
at a set of fc— points in the first Brillouin zone (BZ), and the /i-th Bloch function 
(/»p(fc) is given by 

Mk) = ^Y.^'''''M'--^^^ (4) 
j 

where iV — cxd is the total number of unit cells in the system, (/'^(r — Rj) is 
the atomic orbital (AO) {pz orbital mentioned in section [2T| located in the j*'' 
unit cell defined by the lattice vector Rj . The definition of the Bloch orbital 
corresponding to the spin down electrons will be identical, with a replaced by 
/3 in Eq. [31 Because, in the P-P-P model, the basis functions 0^(r — Rj) are 
assumed to form an orthonormal set, the UHF equation for up-spin orbitals can 
be written in the matrix form as 

F'^ik)C^,{k)^eZik)C:ik) (5) 

where, for a given k value, F°'{k) represents the Fock matrix for the up-spin 
electrons, C"(fc) represents the corresponding c\rn {k) coefficents, arranged in 
form of a column vector, and e"(A:) denotes the band eigenvalue. The Fock 
operator for electrons of up spin is given by 

F"(fc) = h{k) + {r{k) + Jl^{k) - K'^ik)) (6) 

where h{k) denotes the Fourier transform of the one-electron parts of the P-P-P 
Hamiltonian (c/. Eq. [ij, J°'{k)/K'^{k) are the Coulomb/exchange integrals 
for the up spin electrons, obtained by Fourier transforming their real-space 
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counterparts 

= E e'^'^'h^ARj), (7) 

oo 

J^uik) = E e^'^'^^M.l^.), (8) 

oo 

KM = E e'^'^'K^iRj)- (9) 

j=-oc 

Above is the one-electron part of the P-P-P Hamiltonian, and 

J^viRj) / ^p.viRj) denote the Coulomb/exchange mtegrals m the real space, 
defined as 



mm oo 



^M^(^.)-EE E E + 

(10) 



^12 

cr— 1 A— 1 k— — oo 1— — 00 



and 

m Tn oo 



I I V lit- '^A^ ^ 

^".(^.) = EE E ^"A(^fc) E + 

i7— 1 A— 1 fc— — oo — oo 

(11) 

where the expression for a general two-electron repulsion integrals is 
Wi)a{Rj)\^^\u{Rk)\{Ri)) = j j ^^,{vi~R,)dp,{Yi^Rk)r 



-1 
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X c^„{v2 - Rj)cj>x{r2 - Ri)d^ridh2, (12) 
and the density matrix for the up-spin electrons, D'j^^{Rj), is given by 

-, . "a 

DU^^) = X / E C^:(fc)C."„(fc)e^'^-«^dfc, (13) 

n=l 

where the integral over k is performed over the one-dimensional (ID) BZ of 
length A, and ria denotes the number of up-spin electrons per unit cell. Above 
we have given the explicit expressions of various quantities for their up-spin 
components only, because the expressions for the down-spin components can be 
obtained simply by interchanging a and /?. The total energy per unit cell of a 
given system is computed using the real-space expression 

Eceii = E E{^A--(^J>M^(^j) + lD'^^iR^){j;:,{R,) - K^^iR,)) 

+ \K'^{R3){JUR3) - KM) + Dl,{Rj)J<^,{R,)}, (14) 
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where D^i,{Rj) = D'^^j{Rj) + D^^{Rj), denotes the total density matrix of 
the system. The expressions of the two-electron integrals appearing in Eqs. 
[TUl and [TT] are of the most general type, however, in case of the P-P-P model 
(c/. Eq. [T]) only density-density type of e-e repulsion terms are included, which 
implies 

{fi{Ri)a{Rj)\-^\i'{Rk)XiRi)) ^ 5^^5ax5R^RjR^R^V^(^o)x{R,-R,) (15) 

where V^(o)A(fl'j implies that expression is calculated using Eq. [5J assuming 
that the /i-th basis function is located in the reference unit cell while the A-tli 
basis function is in the cell with location Rj — Ri. After the simplification of 
Eq. [ini evaluation of J^^{Rj) and K^^{Rj) becomes quite easy: (a) in Eq. [TU] 
only an infinite lattice sum over Ri needs to be performed, which is done by 
including a large number of terms, and (b) in Eq. [TT]both the sums for Rk and 
Ri reduce to one term each. The convergence of our calculations with respect 
to these lattice sums was tested extensively. 

The UHF equations of the system, leading to the band structure (ei^\k) /e"n \k)), 
and the correpsponding Bloch orbitals, are obtained by solving Eql5l and its /3- 
spin counterpart, by iterative diagonalization technique at a set of fc-points, 
until the total energy per cell of the system (c/. Eq. converges. During the 
self-consistent HE iterations, the integration over the BZ is performed using the 
Gauss-Legendre quadrature technique as suggested by Andre et (author?) 
(H^ l ; with the additional flexibility that the number of points used for the quadra- 
ture can be chosen by the user. 

In order to perform calculations in the presence of a static external electric 
field to simulate the gate bias, one can solve the HE equations using a modi- 
fied Fock operator under the electric dipole approximation by introducing the 
corresponding term containing the uniform electric field E. The modified Eock 
operator of the system is then given by 

Ks^eM = ^" - /^-E = + |e|E.r , (16) 

where F°' is the unperturbed Eock operator for the up-spin electrons in the 
absence of the electric field, e represents the electronic charge, /x = — er, is the 
dipole operator, and r is the position operator for which the usual diagonal 
representation is employed. 

2.3. Density of states 

The density of states (DOS) is obtained using the well-known expression 

p(e) = C^y'e-('-''W)'/2^'dfc (17) 

i 

where e is energy at which DOS is computed, ei {k) is the energy of i-th orbital 
at a given k point, 7 is the broadening parameter, and C includes the rest of the 
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constants. The integration over k is performed over the ID BZ, and summation 
over i includes all the Bloch orbitals of the system. In our calculations we set 
C = 1 to obtain the DOS in the arbitrary units. 



2.4- Theory of optical absorption 

The optical absorption spectrum of the incident radiation polarized in x or 
y direction is computed in the form of the corresponding components of the 
imaginary part of the dielectric constant tensor, i.e., eii{io), using the standard 
formula 

e (co)-C^ f \{c{k)\pMk))\' .18. 

'"^"^ -""i-J {(£;c.(fc)-M^ + 7^}i^cl(fc)*' ^''^ 

where pi denotes the momentum operator in the i-th Cartesian direction, lu 
represents the angular frequency of the incident radiation, E(.v{k) = tc{k) — 
e„(fc), with ec(fc) {e.v{k)) being the conduction band (valence band) eigenvalues 
of the Fock matrix, 7 is the line width, while C includes rest of the constants. 
Assuming that the valence band eigen state \v{k)) is expressed as (c/. Eq. [31 
ignoring the spin orientation) 

h(fc)) =^C^,(fc)|x^(fc)), 

with a similar expression for the conduction band eigen states \c{k)). The 
momentum matrix elements {c(k)\pi\v{k)) needed to compute tuiuj), for a ID 
periodic system, can be calculated using the formula, (author?) [55| 



{c{k)\pMk)) ^ki ^T.^*-cik)C,,ik)-^H,^{k) 

fl.V 



. .mo(ee(fc)-e„(fe)) ^^.^^^^^^^^^^^^,,^ (^g^ 



where Sn implies that the term is nonzero when i denotes the periodicity di- 
rection (a; direction), mg is the free-electron mass, -^H^^^^k) represents the 
derivative of the Hamiltonian (Fock matrix, in the present case) with respect 
to fc, ci{iy denotes the matrix elements of the i-th component of the position 
operator d defined with respect to the reference unit cell, and accounts for 
the so-called intra-atomic contribution(author?) [H^. Note that Eq. [T9l can 
also be used to compute the matrix element {c{k)\py\v(k)) needed to calculate 
the absorption spectrum for the y-polarized light for GNRs (which are periodic 
only in the x direction), by setting the first term on its right hand side to zero, 
and retaining the contribution only of the second term. In these calculations, 
-^H^^{k) was computed numerically, while the usual diagonal representation 
was employed for the d operator. Furthermore, we set C = 1 in all the cases to 
obtain the absorption spectra in arbitrary units. 
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3. Description of the Program 

Theoretical formalism described in the preceding section has been imple- 
mented numerically in a computer program called "ppp_bulk.x" using the For- 
tran 90 (F90) programming language. Advanced features of F90, such as the 
dynamic memory allocation, modules, etc,, have been used to improve the ef- 
ficiency of the code. Architecture independent numerical precision is used for 
portability of the code across different computer architectures. This program 
is capable of doing both the tight binding as well as P-P-P model calculations. 
Using a small number of input parameters such as the positions of the atoms in 
the unit cell, lattice translation vector, hopping and P-P-P Coulomb parameters 
etc., ppp_bulk.x determines the band structure, density of states, joint density 
of states, optical and electro absorption spectrum for ID periodic 7r-conjugated 
systems. As mentioned earlier, the lattice sums are performed in real space 
by including a large number of unit cells, and integration along the BZ was 
performed using the Gauss-Legendre quadrature approach(author?) [HJl. The 
convergence of SCF iterations is slow for systems with large number of electrons, 
therefore, we have also implemented the method of damping to speed up the 
convergence. 

Our computer code consists of the main program, and various subroutines. 
Optionally, the user can link to the LAPACK/BLAS libraries, whose diagonal- 
ization routines can be used by our program. In the following we briefly describe 
the main program, as well as each subroutine/function. 

3.1. Module MTYPES 

In this module the precision of REAL variables used throughout the pro- 
gram is defined. This will facilitate machine independent precision for REAL 
variables. 

3.2. Module MCOMMONDATA 

Global data shared by several routines in the program is defined in this 
module. 

3.3. Main program PPP_ BULK 

This is the main program of our package for performing electronic structure 
calculations within a semiempirical TB formalism for ID periodic systems. It 
calls other subroutines to accomplish various tasks. 

3.4. Subroutine INPUT 

This routine reads input data such as which Hamiltonian to use, its parametriza- 
tion, total number of atoms in the unit cell, their Cartesian coordinates, ID 
lattice constant, and number of /c-points (71^) for BZ sampling etc. Besides, this 
subprogram reads the options to perform various types of calculations such as 
tight binding, RHF, UHF etc. This subroutine also reads the components of ex- 
ternal electric field in the units of V/A, if calculations for a gated configuration 
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need to be performed. Alternatively, by calling other routines, one can generate 
coordinates of some important structural units such as AGNRs and ZGNRs of 
different widths to facilitate an easy realization of the GNR under consideration. 
In addition, various arrays are allocated dynamically and deallocated after the 
desired task is completed. 

3.5. Subroutine GNR_RATOM 

This routine generates the coordinates of various atoms in the unit cell of 
the given GNR, based upon the user specified data consisting of the type of 
GNR (AGNR or ZGNR), its width, and the nearest-neighbor bond length. It 
also computes the lattice translation vector. 

5.6. Subroutine ERROR 

This routine writes out a fatal error message to a user specified logical unit, 
and stops the execution of the code. 

3.7. Subroutine SORT 

This routine, adapted from Ref. (author?) [H^, generates an array in which 
distances between different pairs of atoms in the system are stored in the as- 
cending order. 

3.8. Subroutine GET_NN 

This routine uses the array generated by subroutine SORT to computes the 
distance between pair of atoms which are nearest neighbors (NN), second NN, 
third NN and so on depending the number of unique hoppings defined. 

3.9. Subroutine get_NTij 

This subroutine computes the total number of hopping matrix elements con- 
necting various atomic sites in the system. 

3.10. Subroutine HOPPING 

This routine generates the hopping matrix elements connecting various sites. 
Hopping matrix elements are assigned for a pair of atoms depending on the 
distance between them. The unique hopping matrix elements are defined by 
user in the input file starting from NN atoms, followed by the second NN atoms, 
and so on. 

3.11. Subroutine IJPK 

This subroutine is used to pack the row index i and column index j of 
an element of the upper triangle of a real symmetric matrix into an integer 
corresponding to its location in a ID array. 
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3.12. Subroutine INVIJPK 

The task of this subroutme is just the reverse of the subroutine IJPK, i.e., 
it is used to unpack the integers i (row index) and j (column index) from the 
location of the corresponding matrix element of an upper-triangular symmetric 
matrix packed in a ID array. 

3.13. Subroutine CHEKEDGE 

This subroutine is used to find whether a given site lies on the edge or in the 
interior of the system. It calculates the for the given site, and if NN =2 
then it is regarded as an edge site, and if NN = 3 then it is a site in the interior 
of the system. 

3.14. Subroutine FILOPN 

This subroutine is meant for opening a file which may either be new or old. 

3.15. Subroutine FILCLS 

This subroutine is meant for closing an already open file. 

3.16. Function DOTPD 

This function calculates the dot product between two given vectors. 

3.17. Subroutine GETNA 

Assuming that a given lattice vector (r) is in the form r = na and finds the 
integer n, where a is the primitive vector of the lattice. 

3.18. Subroutine PPP_ PARA 

This subroutine generates parameters associated with the P-P-P Hamilto- 
nian, as per the user choice. 

3.19. Subroutine PRINTR 

This subroutine prints the coordinates of the unit cell in an output file 
called 'unitcell.xsf, which can viewed using the visualization packages such as 
xcrysden(author?) [s^l. Additionally it also generates an output file called 
'system. xsf in which coordinates of atoms in several unit cells are printed for 
viewing the periodic system under consideration. 

3.20. Subroutine CELL_DRV 

This is the driver routine for generating the unit cell related data by calling 
another subroutine named 'cell_ld'. 

3.21. Subroutine GELL_ ID 

Starting with the primitive cell and the lattice vector related data, it gen- 
erates the coordinates of all the cells and atoms included in the calculations. 
This data is useful when lattice sums, to account for the long-range Coulomb 
interactions, are performed. 
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3.22. Subroutine MATEL_R 

This is the master routine meant for generating the one- and two-electron 
matrix elements in the real space, with, or without, the external electric field. 
This is done based upon the data specified by the user in the input file such as 
the Hamiltonian under consideration. Coulomb parameters to be used (if any) , 
hopping matrix elements connecting various sites, etc. 

3.23. Subroutine READ _Ri 

This routine reads the coordinates of atoms in different unit cells, generated 
by the routine cell_ld, for performing real space lattice sums. 

3.24-. Subroutine Makntpq 

This routine arranges the unit cells in an order required for performing real- 
space lattice sums. The reference unit cell is numbered '0', the unit cells to the 
right of reference unit cells are identified labeled with positive integers, and the 
unit cells to its left are labeled with negative integers. 

3.25. Subroutine Maktmat 

This routine stores the hopping matrix elements in a translationally invariant 
format . In this format each hopping element is stored as tij{ipq), where orbital 
i is assumed to be in the cell at location ipq, while orbital j is in the reference 
ceh. 

3.26. Subroutine NUCNUC 

This subroutine computes real space nucleus-nucleus repulsion term for the 
P-P-P model. 

3.27. Subroutine ELE_NUC 

This subroutine computes real space electron-nucleus repulsion term for the 
P-P-P model. 

3.28. Subroutine COULOMB 

This subroutine computes the long-range Coulomb part of the e-e repulsion 
term J^^{Rj) {cf. Eq. for the P-P-P model. 

3.29. Subroutine EXCHANGE 

This subroutine computes the long-range exchange part of the e-e repulsion 
term Kf^i,{Rj) {cf. Eq. [TT|) corresponding to the P-P-P Hamiltonian. 

3.30. Subroutine gen_Kmesh 

This subroutine creates fc-points in the positive part of the first BZ of a ID 
system, between the limits and tt. Various fc-points are non-equidistant, and 
chosen in accordance with the Gauss-Legendre quadrature. The BZ integration 
is performed for the dimensionless variable fca, a being the lattice constant. 



13 



3.31. Subroutine GAULEG 

This subroutine generates the roots and weights needed for the Gauss- 
Legendre quadrature meant for BZ integration. 

3.32. Subroutine SCF_ DRV 

This is the driver routine for carrying out the iterative self-consistent field 
(SCF) calculations. 

3.33. Module MSCF_ VAR 

This module defines the variables common to routines that perform SCF 
iterations. In addition, subroutines for allocation and deallocation of several 
arrays are also defined in this module. 

3.31 Subroutine SCF _RHF 

This subroutine solves the RHF equations for the system under considera- 
tion in a self-consistent manner, using the iterative diagonalization procedure 
at each fc-point, and returns the canonical SCF orbitals, their eigenvalues, and 
the total energy per unit cell of the system. The arrays which are needed during 
the calculations are allocated before the calculations begins, and are deallocated 
upon its completion. Before the first iteration, Hiickel model Hamiltonian is di- 
agonalized to obtain a set of starting orbitals. Subsequently, the Fock matrix 
corresponding to those orbitals is constructed, and diagonalized. The process 
is repeated until the self-consistency is achieved. Depending upon the choice of 
the user, the eigenvalues and eigen vectors can be obtained using either the ZH- 
PEV subroutine from the LAPACK/BLAS library, or the inbuilt HOUSEH_C 
subroutine(author?) [H^, based upon the Householder diagonalization scheme. 

3.35. Subroutine SCF_ UHF 

This subroutine is exactly the same in its logic and structure as the previously 
described SCF_RHF, except that the task of this routine is to solve the UHF 
equation for the system under consideration. Different Fock matrices for the 
a and the /3 spin are constructed and diagonalized in each iteration, until the 
self-consistency is achieved. Eigenvalues and eigen vectors are computed using 
the routine ZHPEV/HOUSEH_C. The iterations are stopped once the total 
UHF energy of the system converges to within a user defined threshold. 

3.36. Subroutine COMMUL 

This subroutine computes the product of two complex numbers, and returns 
the real and imaginary parts of the product, separately. 

3.37. Subroutine GRADFK 

This subroutine finds the fc-space derivative of Fock matrix i^^g^)-, using 
the central difference formula. Recall that is needed to evaluate the 

momentum matrix elements (c/. Eq. [19]) required for computing the linear 
optical absorption spectrum. 
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3.38. Subroutine BAND 

This subroutine computes the band structure of the system by diagonahz- 
ing the converged Fock matrix at a large number of uniformly-spaced set of k 
points. For the purpose, the Fock matrix at those k points is obtained by Fourier 
transforming converged Fock matrix in the real space. The Fock matrix in the 
fc-space is diagonalized using the subroutine ZHPEV/HOUSEH_C to obtain its 
eigenvalues and eigen vectors. 

3.39. Subroutine PRINT _Ek 

In this routine, the eigenvalues of the Fock matrix obtained in the subrou- 
tine BAND are written in various ASCII files. If the calculation was an RHF 
one, then the data is written in the file named 'bands.dat', while for the UHF 
calculations the up-spin eigenvalues are written in the file 'bandsUP.dat', and 
the down-spin ones in the file 'bandsDOWN.dat' . These files can be used 
for plotting the band structure using the standard graphics packages such as 
xmgrace(author?) [5§| or gnuplot (author?) (60| . 

3.^0. Subroutine DOS 

This subroutine computes the density of states (DOS) and writes it in an 
ASCII file named 'dos.dat'. The input to the routine consists of the energy 
windows, and the broadening parameter (c/. Eq. I17p . 

3.U. Subroutine OPTICS 

This is the master routine meant for evaluating the optical absorption spec- 
trum. Several other subroutines are called in this routine to perform specific 
tasks. The range of frequencies over which the spectrum is to be computed, 
along with the line width, are read from the input file. The calculated spectrum 
is written in output files called 'sigma_x.dat' for incident radiation polarized 
along X-axis, and 'sigma_y.dat', for incident radiation polarized along y-axis, 
respectively. Furthermore, the joint-density of stats is written in the file named 
'jdos.dat'. The data in all these files can be visualized using packages xm- 
grace(author?) [H^ or gnuplot (author?) [gOl- In addition, several output files 
named 'sigma_mn.dat', where band-specific optical absorption spectra, due to 
transition from the m*'* valence band to the 7i*''conduction band, are also gen- 
erated. 

3.42. Subroutine SINCOS 

This routine computes coa{kRj) and sin(fci?j) for various values of k and Rj, 
needed for performing the Fourier transforms of various quantities. 

3.43. Subroutine makGmat 

This routine constructs the electron repulsion part of Fock matrix for the 
RHF calculations performed within the P-P-P model Hamiltonian. 
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3.44- Subroutine makGmat_ UHF 

This subroutine is analogous to the routine makGmat, the only difference 
being that it constructs the repulsion matrix separately for the up- and the 
down-spin electrons, and is called when UHF (as against RHF) calculations are 
needed. 



3.45. Subroutine FOUTRA 

This subroutine performs a Fourier transforms on a given complex operator, 
either from the k-space to the r-space, or vice versa, depending upon the value 
of the input variable 'iflag'. Because the systems under consideration are ID, 
the BZ integration is performed only in the positive part of the BZ, as the 
contribution of the negative part of the BZ is same as that from the positive 
part. 

3.46. Subroutine HO USEH_C 
This subroutine, which originally belongs to the EISPACK library(author?) 



58fl , is used for diagonalizing the Fock matrix at different k points to obtain its 
eigenvalues and eigenvectors. For the purpose, the Householder diagonalization 
approach is utilized(author?) [58i |. 

3.47. Subroutine DENK_ RHF 

This subroutine constructs the density matrix in the k space, using the 
Bloch orbitals for the closed-shell systems, assuming that the orbitals are doubly 
occupied. 

3.48. Subroutine DENK_ UHF 

This subroutine performs the same function for UHF calculations, which the 
routine DENK_RHF performs for the RHF calculations. It generates different 
density matrices for the orbitals with a and j3 spins, and also obtains the total 
density matrix by adding them. 

3.49. Subroutine SYMMAT 

This subroutine multiplies the off-diagonal Fock matrix elements of a real- 
symmetric matrix by a factor of 2 to enforce the upper-triangular nature of the 
matrix. 



3.50. Subroutine xMMEcv 

Using ^^jr^ computed in the subroutine 'GRADFK', this subroutine com- 
putes the momentum matrix element {c{k)\px\v{k)) defined in Ea ll9l for incident 
radiation polarized along the x-direction. 

3.51. Subroutine yMMEcv 

This subroutine computes the momentum matrix element {c{k)\py\v{k)) de- 
fined in Eg 1 191 for the incident radiation polarized along the y-direction. 
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3.52. Function LINESHAPE 

This function computes the Lorentzian hue shape, inputs are hne width, fre- 
quency of the incident radiation, and the energy difference between the conduction- 
and valence-band eigenvalues, at a given fc-point. 

3.53. Subroutine SIGMA_ABSORB 

This routine evaluates the optical absorption spectrum defined in Eq ll8l 

3.54. Subroutine printTimeDate 

This subroutines uses the Fortran 90 intrinsic routine DATE_AND_TIME 
to print date and time. 

3.55. Subroutine bloch_ out 

This subroutine generates an output file named 'bloch.dat', in which the 
eigenvectors are written in the binary format. 

4. Installation, input files, output files 

We believe that the installation and execution of the program, as well as 
preparation of suitable input files is fairly straightforward. Therefore, we will not 
discuss these topics in detail here. Instead, we refer the reader to the README 
file for the details related to the installation and execution of the program. 
Additionally, the file 'mcinual . pdf ' explains how to prepare a sample input file. 
Several sample input and output files corresponding to various example runs 
are also provided with this package. 

5. Results and Discussions 

In this section we demonstrate the abilities of our code by presenting results 
for some ID 7r-conjugated systems such as t-PA, PPP, AGNR-11, AGNR-14, 
ZGNR-8 (ZGNR-TV^, denotes a ZGNR with Nz zigzag fines across the width), 
and ZGNR-10. Unless otherwise specified, all the calculations are performed 
using the screened parameters with U = 8.0 eV and Kij — 2.0 {i ^ j) and 
Ki^i = 1. However, later on we examine the influence of Coulomb parameters on 
the band structure of GNRs. In order to benchmark our code, first we use it to 
compute the total energy per unit cell {Eceii ) for t-PA and the polymer PPP, and 
compare it with the results obtained using our earlier P-P-P code meant for finite 
systems (author?) [i^l, with the increasing system size. We also analyze the 
convergence of the total energy/cell with respect to the total number of k points 
(n^) used for the BZ integration, number of cells considered for summation of 
exchange integral {uexc, cf. Eq. [TT|) . by means of calculations on AGNR-6. 
Furthermore, we present our results for the band structure, density of states, 
linear optical absorption spectrum of some of the systems mentioned above. The 
calculations and analysis of the electric filed driven half-metallicity of ZGNR- 14 
and the electro-absorption(EA) spectrum of ZGNR-8 are also presented. 
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Table 1: Variation of total energy (in eV) per unit cell {E^^u) of t-PA and PPP with the 
number of unit cells (A^) obtained using our earlier P-P-P code for finite systems(author?) 
[ZtII . compared with the result for the infinite polymer obtained by the present code (last 
column). 



System 


TV = 5 


iV = 10 


iV = 50 


N = 100 


Infinite Polymer (This work) 


t-PA 


-3.20 


-3.30 


-3.38 


-3.39 


-3.40 


PPP 


-11.66 


-11.73 


-11.79 


-11.79 


-11.81 



5.1. Total energy 

The variation of total energy per unit cell, Eceii for finite fragments of t-PA 
and PPP with the increasing number of cells {N) obtained using our earlier P- 
P-P code for the finite systems (author?) [Z^l is tabulated in Table [TJ while the 
Eceii for the same systems in the infinite polymer limit, obtained by the present 
code ppp_bulk.x, is presented in the last column of the table. The t-PA was 
considered in its dimerized configuration, with the single (double) bond length 
of 1.45A(1.35 A), and the corresponding hopping matrix element to be 2.232 eV 
(2.568 eV). In case of the polymer PPP, the intra-phenyl C-C bond length was 
taken to be 1.4 A with the corresponding hopping of 2.46 eV, while the length of 
the single bond connecting neighboring phenyl rings was assumed to be 1.54 A, 
along with the associated hopping of 2.23 eV. From table[T]it is evident that, the 
there is excellent agreement between the two sets of calculations, which gives 
us confidence about the essential correctness of our code. 

Next, we examine the convergence of Eceii for AGNR-6 with the increasing 
values of k points (n^) used for BZ integration (Fig. [Ija)) and with the increas- 
ing values of the number of unit cells in the exchange sum, Uexc (Fig- IHb)). 
The nearest-neighbor (NN) hopping was chosen to be i = 2.7 eV. From Fig. [l]it 
is obvious that: (a) Eceii converges at about Uk = 25, and remains insensitive to 
further increase in the value rife, and (b) Eceii converges rapidly with Uexe and 
convergence is achieved for riexe = 10. In both cases we have considered about 
10000 cells to evaluate the Coulomb part of two electron integrals (c./ Eg fTO]) . 
For rest of the calculations we have chosen the values : = 50 and Uexe ~ 24. 

5.2. Band structure and the Density of States 

Our code can also be used to perform band structure calculations, along with 
the associated DOS using the TB model, as well as the P-P-P model, employing 
both the RHF and the UHF methods. In what follows below, we present the 
band structure and DOS for a few of the systems, computed using our code. 

5.2.1. Trans-polyacetylene 

The band structure of the dimerized i-PA obtained by the P-P-P-RHF 
method is presented in Fig. [2] (a), with the Fermi energy {Ep) set to 0. The 
band gap occurs at k ~ ir/a, and is obtained to be 2.30 eV, a value in excellent 
agreement with the the HOMO-LUMO gap of 2.31 eV obtained for a finite frag- 
ment of t-PA consisting of 100 unit cells, computed using our earlier molecular 
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Figure 1: Variation of total energy per unit cell of AGNR-6 with respect to: (a) n^, and (b) 
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Figure 2: (a) Band structure, and (b) DOS of dimerized t-PA obtained using the P-P-P model, 
and the RHF method. 



P-P-P code(author?) [13 . The vanishing DOS (Fig. ^h)) around the Fermi 
energy characterizes the band gap, while the DOS peaks near E = ±1.15 eV 
and E = ±6.3 eV are due to the van Hove singularities a.t k = n/a and k = 0, 
respectively. 

5.2.2. Graphene Nanoribbons 

In this section we present the results of band structure calculations on both 
AGNRs and ZGNRs. Furthermore, we also examine the electricity-driven half- 
metallicity predicted in ZGNRs. by computing its band structure in the presence 
of a transverse static electric field. However, first we examine the infiuence of 
Coulomb parameters on the electronic structure of GNRs. In all the calculations 
on the GNRs reported in the remainder of this work, carbon-carbon nearest- 
neighbor distance was taken to be 1.42 A, and all the bond angles were assumed 
to be 120°. As far as the hopping matrix elements are concerned, for both 
AGNRs and the ZGNRs, the nearest-neighbor hopping was taken to be 2.7 eV, 
while, only for ZGNRs, a second nearest-neighbor hopping of t' — 0.27 eV was 
also included. 

Influence of Coulomb parameters. In our earlier work on GNRs(author?) (48| . 
we proposed a set of "modified screened Coulomb parameters," with U = 6.0 
eV and k^j- = 2.0 {i =/= j) and Ki^i — 1, however, we did not present a detailed 
analysis of the parameter dependence of our results. In Tables[2]and[3]we present 
the results of our P-P-P-HF calculations of the band gaps of AGNRs and ZGNRs 
of varying widths as a function of the Coulomb parameter J7, while the screening 
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Table 2: Band gaps of different families of AGNR-JVa {Na = 3p, 3p + 1, 3p + 2, (p > 1, is an 
integer) obtained using our P-P-P RHF approach for U = 6.0 eV and U = 8.0 eV, compared 
with the ab initio DFT and GW results of Yang et at (author?) Il3ll . 



VV IQLIl I mil 1 


Energy gaps (eV) for 3p series of AGNRs 


GW(author?) [IS] 


DFT (author?) [13] 


U = 6.0 eV 


U = 8.0 eV 


0.61 


2.72 


1.13 


2.31 


2.65 


0.98 


2.01 


0.68 


1.39 


2.01 


1.35 


1.68 


0.55 


1.17 


1.63 


Width (nm) 


Energy gaps (eV) for 3p+l series of AGNRs 


GW(author?) [13] 


DFT(author?) [13] 


U = 6.0 eV 


U = 8.0 eV 


0.36893 


5.5 


2.50 


3.29 


3.72 


0.73785 


3.83 


1.60 


2.18 


2.50 


1.10678 


2.83 


1.13 


1.64 


1.90 


1.47571 


2.40 


0.93 


1.33 


1.55 


Width (nm) 


Energy gaps (eV) for 3p+2 series of AGNRs 


GW(author?) [13] 


DFT (author?) [13| 


U = 6.0 eV 


U = 8.0 eV 


0.49190 


1.71 


0.47 


0.41 


0.67 


0.86083 


1.18 


0.33 


0.31 


0.50 


1.22976 


0.94 


0.22 


0.24 


0.40 


1.59868 


0.77 


0.19 


0.20 


0.33 



defined by Kij = 2.0 {i ^ j) and Ki^i = 1 is employed. Furthermore, we also 
compare our results to the ab initio DFT and GW results(author?) [13j. 

An inspection of the two tables reveals the following trends: for AGNRs, 
with U = 8.0 eV excellent agreement is obtained between our HE band gaps 
and the ab initio GW band gaps (author?) [l3| of the = 3p family, however, 
for other families both with /7 = 6.0 eV and U = 8.0 eV our HE band gaps are 
smaller than the ab initio GW results but larger than the DFT results. For 
ZGNRs excellent agreement is obtained between our P-P-P-UHE results and 
the GW results(author?) [l3| for U = 4.5 eV, while with all the larger values 
of U, our approach overestimates the results as compared to the GW ones. 
Thus, we conclude that no single value of Coulomb parameter U provides good 
agreement between our HE results and GW results for all classes of GNRs. 
Therefore, in the absence of any reliable experimental data on the GNR band 
gaps, we have decided to use the original screened parameters of Chandross and 
Mazumdar (author?) [H^j, with U = 8.0 eV and k^j- — 2.0 (i ^ j) and = 1. 

AGNR Band Structure. Before we discuss the band structure of AGNRs ob- 
tained using our approach, we examine the variation of their band gaps with re- 
spect to the three aforesaid families corresponding to the widths A^^ = 3p, 3p-|-l, 
and 3p -I- 2, where p (> 0) is an integer. This classification is based upon the 
TB values of the band gaps which exhibit the relation E^p > E^p+'^ > Ep'+'^{= 
0) (author?) Il0|. However, ab initio density- functional theory (DFT) calcula- 
tions(author?) [10| on these ribbons predicted a different relationship E^P'^^ > 
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Table 3: Variation of ZGNR band gaps with the ribbon width, computed using various values 
of screened Coulomb parameters, and the P-P-P-UHF method, compared with the ab initio 
DFT and GW results(author?) 11311 . Band gaps are also presented at the edge of the Brillouin 
Zone (fc = vr/a), which are more or less width independent. All the gaps correspond to the 
spin-polarized edge states of ZGNRs. 



Width (nm) 


Energy gaps (eV) for ZGNRs 


DFT(author?) [13] 


GW(author?) [13] 


U = 4.5 eV 


C/ = 6 eV 


J7 = 8 eV 


1.136 


0.35 


1.36 


1.34 


1.91 


3.04 


1.562 


0.33 


1.23 


1.14 


1.61 


2.64 


1.988 


0.25 


1.04 


1.00 


1.40 


2.35 


2.410 


0.23 


0.95 


0.88 


1.22 


2.13 


Band gap (eV) at the zone boundary (fc = 7r/a) 




DFT(author?) pj GW(author?) [13] 


U = A.5eV C/ = 6 eV 


C/ = 8 eV 




-0.45 -1.95 


-1.94 


-2.84 


-4.40 



^3p > £;3p+2(-_^ ^^^Yi the important resuh that even for Na = 3p+2, AGNRs 
exhibit nonzero energy gaps, due to the fact that the bond lengths involving the 
edge atoms are shorter than those in the interior. When the decrease in the 
bond length is incorporated in the TB approach by increasing the correspond- 
ing hopping, one also obtains finite gaps for Na = 3p + 2 ribbons, although the 
relation E^p > EgP~^^ > EgP~^^ still holds. Here, we intend to investigate as 
to which of these relationships holds when the AGNR band gaps are computed 
using our P-P-P-RHF approach. In Fig. [3] we present the graph depicting the 
variation of the band gaps of AGNRs for all the three families, with respect to 
their width, computed using the screened Coulomb parameters and our P-P-P- 
RHF method. From the figure it is obvious that while the band gaps of 3p and 
3p + I families are fairly close to each other for a given width, yet the relation 
^isp+i > ^3p > 0), in agreement with the ab initio DFT and GW 

results (author?) [io|, is found to hold for a fairly large range of width. 

Next, in Fig. |H we present the band structure and the DOS of AGNR-14, 
belonging to the 3p + 2 family, computed using our P-P-P-RHF method. It is 
evident from the figure that the AGNR-14 is an insulating system with a small 
band gap of about 0.65 eV, located at the point fc = 0. The DOS presented in 
Fig. Hfb), naturally vanishes in the region of the gap. Furthermore, it exhibits 
several symmetrically placed peaks corresponding to the van Hove singularities. 

ZGNR Band Structure. The nature of the ground state of ZGNRs is quite inter- 
esting, with several authors reporting the existence of a magnetically-ordered 
ground state, with oppositely oriented spins localized on the zigzag edges on 
the opposite sides of the ribbons(author?) [21, .5il, 6l|, a result verified also in 
several first principles DFT calculations(author?) |10l Il4|. As reported in our 
earlier work, P-P-P-UHF calculations performed using the present code also pre- 
dict magnetized edge state as the ground state of ZGNRs(author?) [i^. Here 
we elaborate the underlying physics by performing TB, P-P-P-RHF and P-P-P- 
UHF calculations on ZGNR-10 using the present code, the results of which are 
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Figure 3: (Color online) Variation of the band gaps of AGNRs-A'^Ai for A'^^ = 3p, 3p + 1, and 
3p + 2, families, with respect to their width. All calculations were performed with the original 
screened parameters (U = 8.0 eV and Kij = 2.0 {i ^ j) and ^ = 1) (author?) [s^ l in the 
P-P-P model. 
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Figure 4: (a) Band structure and (b) DOS of AGNR-14 obtained by P-P-P-RHF method, 
obtained using the screened Coulomb parameters. 
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Figure 5: Band structures of ZGNR-10 obtained by: (a) TB method (b) P-P-P-RHF method, 
and (c) P-P-P-UHF method. 

presented in Figs. [S]and[ni At the TB level ZGNR-10 (or a ZGNR of any other 
width) is obtained to be gapless as depicted in Fig. [5^, characterized by flat 
bands near Ep, leading to a intense van Hove singularity at Ep {cf. FiglS^). 
This suggests an instability in the system, and, therefore, a possibility of a 
structural distortion through electron-phonon coupling, or a magnetic ordering 
mediated by Coulomb interactions(author?) jilj. But. the band structure of 
ZGNR-10 obtained by the P-P-P-RHF method (FiglSb) is very similar to that 
obtained by the TB method, except for a small band gap of about 0.25 eV, 
which is an artifact of the RHF approach. Thus, the RHF method by its very 
nature predicts a non-magnetic ground state, even though it takes e-e interac- 
tions into account. However, once we perform spin-polarized calculations using 
the UHF approach which is based upon separate mean-fields for the up- and 
the down-spin electrons, we get the ground state exhibiting edge magnetism, a 
significant band gap of 2.35 eV {cf. FigO)), and the total energy/cell of -55.532 
eV which is lower than -55.006 eV for the non-magnetic state obtained by the 
RHF calculations. 

The fact that ground state exhibits edge magnetism is obvious from the spin- 
density plot for the ZGNR-10 obtained from the UHF calculations, presented 
in Fig. [71 As far as the numerical aspects of the UHF method are concerned, 
it is important to note that while performing the UHF calculations the initial 
guesses for solutions of the up- and down-spin Bloch orbitals must be different. 
Because, for ZGNRs of any width, the number of electrons of spin up- and 
down-spin are equal, as a result of which the UHF solutions converge to RHF 
results, unless the initial guess for the two spin components is different. 
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Figure 6: The density of states of ZGNR-10 obtained by: (a) tight binding method, (b) P- 
P-P-RHF method, and (c) P-P-P-UHF method. The Fermi energy Ep is assumed to be at 
E = 0. 




4 8 12 16 20 
Atomic site 

Figure 7: Spin density of ZGNR-10, obtained using the P-P-P-UHF calculations, plotted at 
different atomic sites of the unit cell across the width of ribbon, starting from top. Anti- 
ferromagnetic alignment of spins across the width is obvious. 
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Figure 8: (Color online) Band structure of ZGNR-14 in the presence of external electric filed 
{Ey = 1 V/nm ), obtained by P-P-P-UHF method. The solid-black (dotted-red) lines represent 
the bands of o(/3)-spin electrons. 



In the presence of a lateral electric field, ZGNRs exhibit half-metallic be- 
havior(author?) |T3|, leading to their possible use in spintronics. In Fig. [8] 
we present the band structure of the ZGNR-14 exposed to a field strength of 
2 V/nni. The calculations were performed using the P-P-P-UHF model with 
i7 = 8. In the absence of the field, as discussed above, the bands of the up- (a) 
and down-spin (/3) electrons are degenerate with a band gap of 1.96 eV . The 
degeneracy is lifted in the presence of electric filed, the band gap for electrons 
of spin a is changes to 1.74 eV and the that for electrons of spin /3 changes to 
0.08 eV, indicating the half-metallic nature. The band gap and the critical field 
strength (i?^) to achieve the half-metallicity for a given ZGNR decreases with 
the decrease in the value of the Coulomb parameter \J . For example, in case of 
ZGNR-14, with I] — 4.5 eV, the band gap for electrons of spin /3(a) changes to 
0.06 (0.76) eV, when the ribbon is subjected to a electric filed strength of 1.0 
V/nm, from the gap of 0.79 eV, in the absence of the field. Furthermore, for 
a fixed value of C/, decreases with increasing the width iw) of the ZGNR, 
due to the decrease in the band gap and increase in the potential difference 
{y — wEy) between the two edges of width i(;(author?) [lOj. 

5.3. Linear optical absorption spectrum 

In our earlier work we argued that the polarization characteristics of the 
optical absorption in GNRs is highly dependent on the nature of its edges, 
and, thus, can be used to determine whether a ribbon has zigzag or armchair 
edges (author?) [S]. In Fig. [9] we present the optical absorption spectrum of 
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Figure 9: (Color online) Optical absorption spectrum of AGNR-14 obtained by P-P-P-RHF 
method. The solid- black line represents e^x, while the dotted-red line represents eyy. A line 
width of 0.05 eV was assumed. 

the AGNR-14 obtained using P-P-P-RHF method. If Emn denotes a peak in 
the spectrum due to a transition from m-th valence band (counted from top) 
to the n-th conduction band (counted from bottom), the peak of exx{^) (black 
line) at 0.69 eV is En, at 3.40 eV is S33. The peaks of iyy{oj) (red/dotted line) 
at 2.05 eV correspond to S12 and E21. The noteworthy points is that individual 
peaks in the spectrum correspond to either x- or y-polarized photons, consistent 
with the D2/1 point group of AGNRs. Furthermore, the x— and polarized 
peaks are well separated in energy. 

In Fig. [TUlwe present the optical absorption spectrum of ZGNR-10 obtained 
using the P-P-P-UHF method. The peaks in exxii-^) are located at 2.41 eV 
(Ell), 3.20 eV (E12), and 4.07 (E22), whereas the prominent peak of eyy{ijj) are 
at 2.41 eV (E") and 3.20 eV (E12). Therefore, in spite of the fact that the 
point group of ZGNRs is also D2h; yet unlike AGNRs, most of the prominent 
peaks of ZGNR-10 exhibit mixed polarization characteristics. This is due to 
the fact that the edge-polarized magnetic ground state of ZGNRs no longer 
exhibits I?2h symmetry because of the fact that the reflection symmetry about 
the a;z-plane is broken, thereby leading to mixed polarizations in the optical 
absorption. Thus, by performing optical absorption experiments on oriented 
samples of GNRs, one can predict whether a given ribbon is AGNR or ZGNR 
by probing the polarization characteristics of the absorption peaks. 

5.4^. Electro-absorption spectrum 

Electro-absorption (EA) spectroscopy, which consists of measuring optical 
absorption in the presence of a static external electric field, has been an impor- 
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Figure 10: (Color online) Optical absorption spectrum of ZGNR-10 obtained by the P-P-P- 
UHF method. The solid- black line represents e^x, while the dotted-red line represents €yy. A 
line width of 0.05 eV was assumed. 

tant experimental tool for probing the electronic structure and optical properties 
of conjugated polymers and other materials(author?) [6^. EA spectrum is de- 
fined as the difi^erence of the linear absorption spectra with, and without, an 
external static electric field. In our earlier work we argued that the EA spec- 
trum can be used as a probe of both the electric-field driven half-metallicity, 
as well as the edge magnetism of ZGNRs (author?) [i^. The essential physics 
behind it is that in the presence of a lateral external electric field, the band gap 
of a spin-polarized ZGNR for spin-up electrons is different from those of down 
spins, leading to two split optical transitions across the gap. We illustrate this 
in Fig. [Tl] which contains the EA spectrum of ZGNR-8 in the presence of a 
lateral external electric field of strength 1 V/nm. as well as the linear absorp- 
tion without the field, for its spin-polarized ground state, computed using the 
P-P-P-UHF approach. The tendency towards half metallicity is apparent with 
the presence of two energetically split peaks corresponding to two different Sn 
transitions among up and down-spin electrons. Noteworthy point is that the 
two split peaks across the fundamental gaps will occur only if the ZGNR has 
a edge-magnetized ground state, and not when the system has a non-magnetic 
ground state. Therefore, if these split peaks predicted in our work can be ver- 
ified in the EA spectra of the ZGNRs, it will provide an all-optical probe of 
determining the nature of the edges in GNRs. 
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Figure 11: (Color online) Linear absorption spectrum (solid-black line), and the electro- 
absorption spectrum (dotted-red line) of ZGNR-8 obtained by the P-P-P-UHF method. 

6. Conclusions and future directions 

In this paper we have described our Fortran 90 program which solves the HF 
equations for both the closed- and open-shell ID-periodic 7r-conjugated systems 
using the TB and P-P-P models. The present computer program has been writ- 
ten in Fortran 90 language which allows dynamic allocation of memory, thereby 
freeing the code from artificial limits related to array sizes. To demonstrate the 
capabilities of our code, we presented results of numerous test calculations on 
various systems including organic polymers, as well as GNRs. We presented 
the results of total energy calculations on the polymer t-PA and PPP, while 
the band structure and the density of states of t-PA and various GNRs were 
reported. Furthermore, we also explored the electric-field driven half metallicity 
of ZGNRs, as also the optical absorption spectra of both AGNRs and ZGNRs. 
We also reported calculations on the EA spectrum of ZGNRs, and argued that 
it can be used to determine the nature of edge termination in GNRs so as to 
differentiate between armchair and zigzag type edges. 

Having developed a HF mean-field code for the ID periodic 7r-conjugated 
system, our next aim is to include the electron-correlation effects for them, 
particularly for the band structure, as well as to account for the excitonic effects 
in optical absorption. Work along those directions is underway in our group, 
and results will be communicated in future publications. 
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